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Abstract 

The analytic expression of the time evolution of the Reynolds stress 
anisotropy tensor in all planar homogeneous flows is obtained by exact inte- 
gration of the modeled differential Reynolds stress equations. The procedure 
is based on results of tensor representation theory, is applicable for general 
pressure-strain correlation tensors, and can account for any additional turbu- 
lence anisotropy effects included in the closure. An explicit solution of the 
resulting system of scalar ordinary differential equations is obtained for the 
case of a linear pressure-strain correlation tensor. The properties of this so- 
lution are discussed, and the dynamic behavior of the Reynolds stresses is 
studied, including limit cycles and sensitivity to initial anisot ropies. 

I. INTRODUCTION 

The dynamical behavior of the Reynolds stresses in homogeneous flows is modeled by a 
tensor evolution equation. Previous studies have focused on the fixed points associated with 
the equilibrium states of the Reynolds stresses for several homogeneous flows 1-3 in order 
to assess the stability of higher order models and the ability of these models to reach the 
correct solution points. In recent studies, such fixed points have been obtained analytically 
for all planar homogeneous flows in both inertial and noninertial frames as asymptotic states 
of the evolution of the Reynolds stresses . I. * * 4 In the present paper, the time evolution of the 

Reynolds stress anisotropy tensor is obtained analytically for all planar homogeneous flows. 
The resulting explicit expression for the Reynolds stress anisotropy tensor is quite compact 
and can be expressed as ratios of sums of exponentials in time. 

Such an analytical solution is obtained through a recasting of the tensor equation for 
the Reynolds stress anisotropy into an equivalent set of three scalar ordinary differential 
equations in three scalar invariants by using representation theory. This procedure can be 
applied to the Reynolds stress model equations in which the pressure strain correlation ten- 
sor is modeled in a general way, including quadratic or higher order terms, and additional 
anisotropy effects can be incorporated. The solution of the resulting set of ordinary differen- 
tial equations is obtained for the case of a linear pressure-strain correlation tensor, with no 



additional anisotropy effects included. The present explicit nonequilibrium stress solution 
predicts stress anisotropies that are quite close to the ones given by the modeled Reynolds 
stress anisotropy evolution equation over all times. The differences can be attributed to 
the assumption of a slow variation of the relative strain parameter that had to be made in 
obtaining the explicit expression of the stress anisotropies. All the dynamic features of the 
Reynolds stress evolution are captured by the explicit time solution, including limit cycles. 


II. EVOLUTION OF REYNOLDS STRESS ANISOTROPY 

Consider incompressible, homogeneous turbulent flow, where the velocity and the 
kinematic pressure p are decomposed into the ensemble mean and fluctuating parts: 


Hi = Ui + u' i , p = /) + //. 


( 1 ) 


In homogeneous conditions, the velocity gradients dui/dxj are independent of position. 
These gradients are also assumed to be independent of the time. The Reynolds stress tensor 
Tij = t/'t/' is a solution of the time evolution equation 


Tij = - tv * 7 


duj diii r 

7T U k "X h Q ij ( t-nik j ^~ik T ^niki^jk ) t 

OX k OXk 


( 2 ) 


which is valid in an arbitrary noninertial reference frame that can undergo a rotation with 
angular velocity Q m relative to an inertial frame. In (2), e,j*. is the permutation tensor and 


, ( dll': fhl'j 


2 ,y < dl ^ dU 'l 


dx h dx k 


(3) 


are the pressure-strain correlation and the dissipation rate tensors (where v is the kinematic 
viscosity ) , respectively. 

With the turbulent kinetic energy K = the scalar turbulent dissipation rate 

e = hen. and the Reynolds stress anisotropy tensor 


b = Tij 


i 


13 ~ 2/\ 3^’ 


(4) 


the term ^ is modeled in the commonly used second-order closure models in the general 
form 2 as 

*ij = -e(C? + C{£) b tJ + C 2 I\$ij + C 3 K (bikSkj + S tk b k] - § b mn S nm S tJ ) 

C . 4 A ( b{ k If k j II i k b k j ) C' 4 1\ ( b{ k ( mk j € mi k b k j ) (5) 

f 5^ (b; k b kj jb mn b nm bijj , 

Above, the strain rate S tJ and rotation rate l V t j tensors are defined as 


c 1 ( chh duj 

1 13 2 \dxj + dxi 


W t , 


1 ( du t duj 


2 \dxj dxi ' 


( 6 ) 
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and V = —T,jSij = —'IKbijSij is the turbulence production. The coefficients Cj\ ( '* , and 
C b • ( V> can, in general, be functions of the invariants formed on b tJ , S,j. and If . Equation 
(5) can be shown to be the most general form for <!>,,. For example, the pressure-strain 
model of Speziale, Sarkar and Gatski 2 (SSG) gives the following coefficients: 


C° = 3.4, Cl = f .8, C 2 = 0.36, C 3 = 1.25, C 4 = 0.4, (\ = 4.2. 


The substitution of (5) into (2) yields the following general evolution equation for the 
Reynolds stress anisotropy tensor b tJ . written in matrix form and in nondimensionalized 
variables 

4-b = - — b - aj bS* + S*b - |{bS*}l) + « 2 (bW* - W*b) ( 8 ) 

at* gij \ •> / 

+ ja 4 (b 2 - l -{b 2 }l) - aiS* -L\ 


where 

S* = S /v ^}. W* = w //{&}. L* = L/V^}, (9) 


V = Tyf{&}, C = Ty/-{W*},r = tyf{&}. (10) 

The tensor accounts for noninertial effects 

h ij — i i ij t in ij i (11) 

where c w = ( C 4 — 4 )j{C\ — 2). The following definitions are used for the coefficients: 

r/’i + 1 ) 7 + ^e“- 1 . (i 2 ) 


{J 


«, = - (' 2 ) /2, «2 = (2 - C 4 )/ 2, « 3 = (2 - C 3 )/2, «< = C 5 /2. (13) 

The tensor L can generally contain the additional turbulence anisotropic effects, and the 
scalar coefficients a t may generally be functions of the invariants ?/ and (. In the current 
context , L is taken simply as d to represent the effects of the dissipation rate anisot ropy, 
with the dissipation rate anisotropy tensor defined as 


Equation (8) is equivalent to (2), but must be supplemented with an equation for the tur- 
bulent kinetic energy A 

K = V — e, (13) 


and for closure, an equation for the turbulent dissipation rate 


-2 

t = Cel 7 - P — Cs 2 ~rr. 

A A 


where C f i and C s2 are closure constants. 


(16) 
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III. SOLUTION OF REYNOLDS STRESS EQUATION 


A. Equivalent Scalar Representation 


The tensor relation (8) governing the evolution of the stress anisotropy cannot be manip- 
ulated further because it involves matrix products and their transpose. Even with lineariza- 
tion (a 4 = 0), the terms that factor b cannot all be grouped to allow for the integration of 
the system of ordinary differential equations. The following technique, however, transforms 
the tensor relation into an equivalent system of scalar ordinary differential equations, which 
in turn can be solved. 

With the evolution of the anisotropy tensor b governed by equation (8), the tensor b can 
be assumed to be dependent only on the tensors S * and W*, as well as on scalar quantities 
such as f*, r/, and It can be shown 5 in this case that for two-dimensional flows the exact 
representation for the tensor b is given by 

fbW*S*l / 1 \ 

b = {bS*}S* + 1 ; (S*W* - W*S*) + 6{bS* 2 } ^S* 2 - -I J . (17) 

Equation (17) thus shows that if the three scalar invariants {bS*}, {bW*S*}, and {bS* 2 } 
can be determined independently of (17), then a knowledge of these scalar functions is 
equivalent to knowing b. In addition, the representation (17) can be used to construct (see 
Appendix A) the nonlinear term in (8), 


b 2 - l{b 2 }I = 2{bS* }{bS' 2 }S* + 2 ^ bW ^J bS ^ (S*W* - W*S*) 


: is) 


+ 


< bS ’> 2 - 2 <b ( ^y - HbS- 2 }^ (S‘‘ - JI) , 


which clearly shows the same tensor function representation as in (17), as well as a depen- 
dency on the same three scalar invariants. Independent of the representations shown in (17) 
and (18), equations for the three scalar invariants {bS*}, {bW*S*}, and {bS* 2 } can be 
formed (see Appendix B) from the Reynolds stress anisotropy evolution equation given in 
(8). For simplicity, the following variables are introduced 


B x = {bS*}. B 2 = {bW*S*}. B 3 = {bS* 2 }, (19) 

and the representation in (17) is rewritten as 

b = B x S* - |f(S*W* - W*S*) + 6 B 3 (s* 2 - ll) , (20) 

where 

3 = e iw- 2 } 

?/ 2 { S * 2 } 

The evolution equation (8) is, therefore, equivalent to the system of ordinary differential 
scalar equations in the scalar invariants {bS*}, {bW*S*}, and {bS* 2 } that is obtained as 
shown in Appendix B (equation (B5)), 
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( 21 ) 


B i — (2a ft — — )ft + 2aj ft> — 2q.->,B$ + - — -ft.ft — o \ — L\. 

V V 

ft = - a 2 H 2 B , + (2a B i - -)B 2 + —B 2 B :i - L 2 , 

V V 

ft = 4« 3 ft + (2«ft - -)ft + ^-ft 2 + -^7? 2 - -S 2 - I 

7/ 0 7/ 47/ /V. 7/ 

where the tensor L appears through the invariants 

U = {L*S*}, />2 = {L*W*S*}, = {L*S* 2 }. (22) 

Equation ( 21 ) is a system of three algebraic ordinary differential equations in the three 
unknowns B \ , B 2 , and # 3 , which is quadratic even if a..\ / 0 because <j depends on B\* 

— = —2ai]B\ + ;i , (2d) 

fj 

with o = (2/ /2+ 1 and T = C^j'2 — 1. Note that t he degenerate case of // = 0 is not considered 
because either the absence of mean velocity gradients or the absence of a turbulence field 
would be implied. Of course, no such restrictions apply to so the case of ( = 0 is not 
precluded. 

The dynamic system (21) is subjected to the initial conditions 

ft,o = {b () S*}. ft. 0 = {b«W*S*}, ft. 0 = {b () S' 2 }. ( 21 ) 

where b (J is a given initial anisotropy. 

B. Solution of Reynolds Stress Equation 

The system of scalar ordinary differential equations (21 ) with the representation in 
(20), which is equivalent to the original tensor evolution equation for the Reynolds stress 
anisotropy ( 8 ). has a significant advantage in that it is much more tractable and better suited 
for analysis than the original tensor equation. Any expression for the extra anisotropy tensor 
L can be provided which involves the stress anisotropy tensor to any degree of complexity. 
It then suffices to study the resulting dynamical system (21) to have a complete description 
of the evolution of the Reynolds stress anisotropy tensor ( 8 ). In the case of pressure-strain 
rate models that are only linear in the Reynolds stress anisotropy so that r / 4 = 0 (compare 
(5) with = 0) and for which no additional anisotropies are included L = 0, an explicit 
solution of the system of ordinary differential equations (21) can be obtained. The solu- 
tion procedure for the resulting differential system is not straightforward, and the major 
steps of its derivation are given in Appendix C. The final expression for the Reynolds stress 
anisotropy tensor, which is the solution of the modeled evolution equation ( 8 ) with « 4 — 0 
and L = 0, is rather compact and involves ratios of characteristic functions ty;: 
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(25) 




«3 
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The characteristic functions are the fundamental solutions of a quadratic nonlinear system 
of two ordinary differential equations (see Appendix C) and are related by and 

^3 - ^2, 

¥t(/*) = A' [ELl /'r^ 4 * + V{H + //°)] , 

«2(n= A’E? s )Prt Ar, ‘, ( 2fi ) 

« 3 (n= AEr=l^rA r e Ar '*, 

where 

H = (27) 

//° = 40(02/^2,0 “ ^3^3,0), (28) 

// r = [A* - 2 aXrB w ~ (H + H°)](\ p - A,), (29) 

1 / = (A3 — A2 )/ Ai T ( Ai — A 3 )/A2 + ( A2 — Aj )/ A3, ( 30 ) 


and A' = [(A 2 — Aj )( A 3 — At )( A 3 — A 2 )] — 1 • In (29), the indices p and q are such that e rpq = — 1. 
Finally, the parameters A r are the eigenvalues that are obtained as roots of the following 
third-order characteristic polynomial equation: 

A 3 — —A 2 — (H + 2aa\ ) A + —H — 0. (31) 

V V 

In (25), the relative strain parameter depends on the time 77 = and its evolution is 

governed by an additional equation. However, in the derivation of the explicit solution of the 
system of ordinary differential equations (21), the relative strain parameter 1 7 was assumed 
not to vary in time, 77 ~ 0. (See Appendix C.) 


IV. ANALYSIS OF DYNAMICAL BEHAVIOR 
A. Transient Behavior 

Equations (20) and (25) completely determine the solution of the modeled evolution 
equation for the Reynolds stress anisotropy tensor for all planar homogeneous turbulent 
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flows. Any initial stress anisotropy can he taken into account. The scalars B, of the ex- 
pansion in (20) involve ratios of the characteristic functions which are expressed as the 
sum of three exponential functions. Because the arguments of the exponentials are the same 
for all characteristic functions and are given as the roots of the characteristic polynomial 
(31), the dynamical behavior of the stress anisotropies will essentially he determined by the 
location of these roots in the complex plane. If Ai and A 2 are defined as 

1 ( 1 iP\ 1 ll ( 1 lf 2 \ 

X . -- (ff + 2na, + A,= - ~ (« - ■ TO 

then the discriminant of the third-degree polynomial equation (31) is given by 

A = A, ! + A 2 . (33) 


The discriminant, depends on the parameters // and // only and can he rewritten as 


A = - 


2 1 


ii A i 2 

H— + (a 2 a 2 ' ,n " ■*” 2 


, + 10 cm i II — 2 II ) — T ( H + 2oy/i ) 

v 4 r 


(34; 


Because H is a function of 7 Z (see equation (27)), the value of the discriminant will he 
essentially determined hv 7v and ?/. Figure 1 shows the evolut ion of the discriminant A as a 
function of rj. for different values of the parameter 7 Z. Three cases must he distinguished: 





FIG. 1. Evolution of discriminant of cubic root equation as function of // for different values of 
parameter 7 vh as labeled. 


(a) A < 0 

The three roots of the characteristic polynomial are real, and the characteristic lunctions 
are combinations of real exponentials. Because a property of the roots of a third-degree 
polynomial is that A! + A 2 + A 3 = ii/t) is always positive, at least one root is positive. 



(b) A > 0 


Two roots are complex conjugates, for example, A 2 = A 3 = d + iw and \ x — A. The 
characteristic functions can then be expressed as 

^i(r) = K f \e xt * + e dt *{f n cos u>t* + f X2 sinu A*) + / 13 ] , 


*2(n 

W) 


K' 

A" 


\e M * + e dt *(f 2 i cosutf* + / 22 sinu?/*)| , 
X 2 e xr + e dt *{f 3 i cos u>t* + f 32 sin a;**)] , 


(35) 


where A ' and /q are constants. Because a property of the roots of a third-degree polynomial 
is that A(c/ 2 + u/ 2 ) — and A + 2d = / 3/r / always, A will be positive when II < 0. 

When H > 0, A < 0 and, thus, 2d = /i/r/ — A > 0. Therefore, in this case also, one root will 
have a positive real part. 


(c) A = 0 

The roots are all real, and two of them are equal, for example, A 2 = A 3 = A. In this 
case, /q = 0, and the third root X x has no effect on the characteristic functions. Because 
2A + Aj = jj/r) > 0, at least one root will be positive. 

In summary, at least one root will always have a positive real part; therefore, the charac- 
teristic functions will always grow exponentially. In the case of A > 0, the characteristic 
functions grow with superimposed (damped) oscillations of period T — 2; r/u;. Three distinct 
cases can be identified on Figure 1. For mean flow fields such that H > 0, that is, 

\iz\ < n u 

with 


(e.g., |' R\ < 0.271 for the SSG pressure-strain correlation coefficients), the discriminant A 
is always negative, and the three roots will always be real. For values of H such that 
— 2 qui < H < 0, that is, 



7vi < 1 72. | < 

(e.g., 0.271 < \R\ < 1.231 for the SSG), where 

_ i [TZ 

K 2 = — \ -«3 + Ofl ], 

«2 V 3 

the roots will have a different nature depending on the magnitude of T], and the evolution 
of the stress anisotropy will contain an oscillatory component for sufficiently small values of 
7/. Finally, for values H < — 2aaj, that is. 
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\n\ > ti- 2 . 

the discriminant A is always positive, and the evolution of the stress anisotropy components 
will contain damped oscillations. 

The case of a vanishing root is of particular interest. Because A) A 2 A 3 = —HjJftf is 
always verified, one of the roots is zero (e.g., A 2 ) if either H = 0 or 1 /// = 0. In the case 
H = 0 (and, thus, % = Tvi), A < 0 always, and the nonzero roots are real and given by 
A !„•{ = 8/{2r]) ± \J[ft/(‘2i})] 2 + 2 cui\. In the case 1/r/ = 0, A = — (// + 2<v<ii ) 3 /27 . When 
H > — 2a«i (i.e., 7v < R 2 )« the roots are all real and are given by the same relations as 
for the case H = 0. When H < —'ion, (i.e., Tl > 1Z 2 ), A > 0, and the roots are purely 
imaginary. The characteristic functions T 2 and have a purely oscillatory behavior: while 
is increasing as / ty 2 . The period of the oscillation is T = 2 tt /A . with us 1 = —{11 + 2o«i). 


B. Asymptotic States 


It has been shown in the previous section that the characteristic functions T, are always 
increasing, except when 1/r/ = 0. For 1/r; / 0, when the effect of the initial conditions has 
vanished, the exponential that corresponds to the root with the largest real part becomes 
dominant, and the ratios of the characteristic functions converge to the values 


where 


* 3 (n = j. = 


jMn ] = 1; jMn = j_ 

^ 2 (/*)J 0C , W) a x 


(36) 


X = max R( { A r ). 

r=l,2.3 


The asymptotic values of the coefficients arc given by 




^-(— - Aoo), 

2o r/, x , 

20 7/oc, 

^-d -T-— )- 
b« Aoo */x 


(37 


where ?/, Xl is the equilibrium value achieved by the relative strain parameter. It can be shown 
that for any planar homogeneous flow described by the Reynolds stress model equation (2), 
a unique relationship holds between the equilibrium value lor the production-to-dissipation 
ratio (^), x ,, the equilibrium relative strain parameter 7/, x , and the rotation rate R: 





2g, 


a i 


V 


-i 




lor — 7v|iin < R< Rhrn- 


— = 0. otherwise 

l 'lL 


( 38 ) 
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where g, x , = (a(r)<x> + fi) 1 and 


^-lim = — 

a 2 



+ aaj + 


0a t 



(39) 


The equilibrium value of the production-to-dissipation ratio (j)oc, is determined by the K 
and f evolution equations. For the standard approach, where equations (15) and (16) are 
used, this value is given by 



a 


e2 


— 1 



(40) 


For values of the parameter 7?. outside the range [— the asymptotic value of the 
relative strain parameter is l///<x, = 0, and the representation coefficients given by (37) are 


B? = — A 00 /(2 o), B? = a 2 n 2 /(2a), B? = a 3 /(6o). (41) 


However, as shown before, for values of l/rj = 0, the solution reaches a limit cycle for the 
anisotropy, and no asymptotic state exists. The solution (41) is, therefore, spurious because 
the real behavior of the anisotropy is purely oscillatory in time. 

In a previous study, 6 an expression equivalent to (37) was obtained from a direct analysis 
of the asymptotic state of (8). Written in the present formalism, the asymptotic value for 
the representation coefficient was obtained from the roots of a cubic polynomial in B^, 

4a 2 (B?f - 4 a-(B?f + ( l — - H - 2aa, | B? + = 0, (42) 

v W J n 

which led to the problem of choosing one of the three roots so that the value of B™ was 
retained. This question could not be rigorously answered, and the selection of the proper 
root was done on the basis of continuity arguments. 6 With the present dynamic approach of 
the Reynolds stress equation, the proper choice for the roots in (42) is obvious and is based 
on the limit of a dynamical process. Clearly, the correct root is the one that controls the 
asymptotic behavior of the system (i.e., A^,). In terms of B^\ B™ must be taken as the 
root in (42) that has the lowest real part. 

In general, planar homogeneous flows can be described by the expression' 

~dx~ ~ 2 + ( D ~ u)hi2&3\\ - (43) 

which yields 

V 2 = \{Brf, C 2 = ^ [(w - 2o w fi) t] 2 , (44) 

where D/'2 is the strain rate and ui/2 is the rotation rate of the flow. As shown in Table I, a 
wide class of homogeneous flows, both with and without system rotation, can be described 
in terms of Tv, 
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TABLE I. Characterization of Common Homogeneous Turbulent Flows. 


Flow 

\*/D 1 

\n/D\ 

mi 

II 

Plane shear 

1 

0 

1 

-1.19 

Plane strain 

0 

0 

0 

0.09 

Hyperbolic* 

< 1 

0 

< 1 

> -1.19 

Elliptic* 

> 1 

0 

> 1 

< -1.19 

Rotating plane shear 

1 

0.25, 0.50 

0.125, 1.25 6 

0 .07. -1.91 


a See Leuchter and Benoit ' for a. description of this class of flows. 

^ These' values are dependent on the pressure-strain rate model used (SSC model in this case). 


7v 2 = 




2 


(15) 


As discussed above, the value A = 0 divides the plane (7?.. //) into two regions in which 
the Reynolds stress components have distinctly different behaviors in time. These regions 
are illustrated in Figure 2; the solid lines are determined by the locus ol points ( 7v . ?/ ) such 
that A = 0. For realizability (// > 0), only the positive root of \/ij is plotted, figure 2 
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FIG. 2. Map of time-evolution types for Reynolds stress tensor in (7v, 1///) plane. Boundary 

A = 0, ; locus of equilibrium points 1/?/^, as function of 7v, ; plane shear, Q; plane 

strain, O; rotating plane shear 1^1 = 0.25 and 0.50, A and Vi respectively. 


also shows the locus of asymptotic solutions l/r/ rX , as a function of 7£, as defined by (38). 
(Note the dashed line in Fig. 2.) The symbols correspond to several planar homogeneous 
flows. (See Table 1.) When the standard equation (16) for the dissipation rate f is used, the 
evolution equation for the relative strain parameter // 
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(46) 


= 2t}Bi(C s \ — 1) + (C £ 2 — 1) 

is solved in conjunction with the evolution of the coefficients B{. Because for a given planar 
homogeneous flow the value of the parameter 7Z is fixed, the system will evolve along vertical 
lines in Figure 2. For values of (7Z,rj) situated in region I of Figure 2 , the roots of the 
characteristic polynomial are real, and the Reynolds stress components converge to the 
asymptotic solution as ratios of real exponentials. In terms of dynamical systems, the 
asymptotic solution is a sink. For example, planar strain flows and rotating shear flows with 
il/D = 0.25 will always have an evolution that is characterized by growing exponentials, 
for any initial condition on the anisotropy b 0 or on the relative strain parameter tjq. Points 
in region II have a time evolution with a damped oscillatory character, and the asymptotic 
state is a spiral sink. The rate of damping of the oscillations is proportional to 1 /?/, with no 
damping at all when \/t) — 0 . For example, a shear flow with high rotation (f l/D = 0 . 5 ) 
is such that 7v > 7v 2 , and the stress components will evolve to their asymptotic value with 
damped oscillations. Note that lor the homogeneous shear case < 7v < 7 \ _>) . the two 
types of evolution can be experienced depending on the value of ij. As already mentioned, 
lor values 1Z > 7ci, tn . the asymptotic solution for the relative strain parameter is I / ?/ i>; = 0, 
and the solution is purely oscillatory, i.e.. a limit cycle is reached. 

V. ILLUSTRATIONS 

First, consider a sheared flow (7v = 1 ) in which the turbulent field is subjected to tin* 
following initial conditions: 


//0 — 4.38, &n,o — ^ 12,0 — ^ 22 , u — 0 . 

Figure 3 show's the evolution in time of the stress anisotropies predicted by the differential 
equation ( 8 ) and by the present, explicit time solution. Clearly, the anisotropies given by 
the present explicit solution are almost indistinguishable from those given by the differential 
equation; the difference is attributed to the assumption that dij/dt* ~ 0 , which is used 
in deriving the explicit solution. (See Appendix C.) From the standpoint of a dynamical 
system, it is more interesting to consider the evolution of the system variables // and b tJ in 
the phase plane, as showrn in Figure 4 . 

In the case of an initial anisotropy, for instance, 

'/o — 3.38, b \ i u = — 0.1, f*i 2 ,o = 0.2, 622,0 — 0.2, 

the present explicit nonequilibrium solution leads to stress evolutions that are almost indis- 
tinguishable from those obtained with the differential Reynolds stress equation, as shown in 
Figure 5. Moreover, the explicit nonequilibrium solution is remarkably close to the differen- 
tial stress equation over a wide range of initial values tj 0 for the relative strain parameter, 
as illustrated in Figure 6 , which shows the initial value of ?/ varying from 1 to 100, with 
isotropic initial conditions ( b , , 0 = 0 ). 

For values of the parameter 7Z outside the range [ — Tvij ril , 7v|j In ] . the asymptotic value 
for l/?/ is shown to be 0 (see equation (38)). and the solution reaches a limit cycle for the 
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0.10 / 



-0.20 l--* — ■— * * 1 ‘ * ‘ 1 * * * — 1 

0 5 10 15 20 

t* 

FIG. 3. Time evolution of stress anisotropies for homogeneous shear case. Initial conditions are 

//o = 3.38; ft] i.o = fti 2 ,o = 622,0 = 0. Present nonequilibrium solution. ; differential Reynolds 

stress equation* . 



U 

FIG. 4 . Phase plane evolution of stress anisotropies for homogeneous shear case. Initial con- 
ditions are 770 = 3.38;6n.o = 612.0 = 622,0 = 0 . Present nonequilibrium solution, : differentia] 

Reynolds stress equation, - • asymptotic solution, o. 

anisotropy. Figure 7 shows the time evolution of the stress anisotropy components in the 
case of a rotation-dominated flow for which u)jD — 2 (and, thus. TZ — 2 ), which is well 
outside the range [— lZ\\ m , 7vii, n ]. The discriminant A is, therefore, always positive. (See 
Figure 2.) The initial stress field is taken to be isotropic, and the initial value of the relative 
strain parameter is arbitrarily set to a high value (7/0 = 100 ) in order to show clearly the 
characteristic oscillation of the dynamic system. From its initial bounded value, the relative 
strain parameter ;/ grows unboundedly (so that 1 / 7/ — > 0 ) with superimposed oscillations, as 
illustrated in Figure 8 . 

In Figures 7 and 8 , clearly the present nonequilibrium explicit solution is extremely 
accurate in capturing the initial phase of the evolution of the anisotropy. The period of the 
oscillations is also captured well by the present nonequilibrium solution. When 1 /// = 0 , the 
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T1 

FIG. 5. Phase plane evolution of stress anisotropies for homogeneous shear case. Initial condi- 
tions are // o = 3.38; 6n,o = —0.1; 6j2,o = 0.2; 622,0 — 0.2. Present nonequilibriurn solution, ; 

differential Reynolds stress equation, ■ • asymptotic solution, 0. 



FIG. 6. Phase plane evolution of stress anisotropies for homogeneous shear case. Initial con- 
ditions are 61^0 = 612,0 — 622,0 — 0, and different values for tjq are used, as labeled. Present 
nonequilibrium solution, ; differential Reynolds stress equation, ; asymptotic solution, o. 

solution is purely oscillatory, as discussed before. Because 7/ starts from a bounded value, the 
oscillations start with a damping component, and their amplitude first decreases in time. 
Although the frequency of the oscillations is captured well, the amplitude clearly is not 
correctly represented by the present nonequilibriurn solution for larger times. This finding is 
attributed to the hypothesis dij/dt* ~ 0 that is used in the solution procedure. As the initial 
condition r/o takes lower values, the initial damping of the oscillations is stronger, so that as 
the limit cycle is approached all oscillations may be nearly killed for the differential stress 
evolution; whereas the oscillations for the present explicit solution have not been damped 
fast enough, as illustrated in Figure 9, where the initial // value is set to a low value (;/ 0 = 2). 
As rj grows in time, the oscillations of the nonequilibriurn solution are not damped at the 
correct rate, and the remaining long-term amplitude of the oscillations is not correct. 
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0 20 40 60 80 100 

t* 

FICi. 7. Time evolution of stress anisotropies for rotation-dominated flow ( / D — 2). The 

initial conditions are t ? 0 = 100; b n , 0 = 612,0 = l>22, o = 0. Present nonequilibrium solution. : 

differential Reynolds stress equation, . 



FICI. 8. Phase plane evolution of stress anisotropies for rotation-dominated flow ( / 1) = 2). 

Initial conditions are 7/ 0 = 100; 6 n ,u = 6 12 .o = 622,0 = 0. Present, nonequilibrium solution. : 

differential Reynolds stress equation, . 

VI. CONCLUSIONS 

A general procedure has been developed that allows for the investigation of t he tilin' 
evolution of the Reynolds stress anisotropy components in all planar homogeneous turbulent, 
flows. Fhe procedure takes the evolution equation for the Reynolds stress anisotropy tensor 
and replaces it with an equivalent system of scalar ordinary differential equations. This 
equivalent system can then be used for assessing the dynamical behavior ol a variety of 
turbulence closure models. This includes pressure-strain rate models which are quadratic 
(or higher) in the anisotropy tensor and in which other anisotropic effects, such as dissipation 
rate anisotropy, can be taken into account. For the case of linear pressure-strain rate models, 
the system of ordinary differential equations can be analytically integrated when the relative 
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FIG. 9. Phase plane evolution of stress anisotropies for rotation-dominated flow (lj/D = 2). 

Initial conditions are 7/0 = 2; 6n,o = ^ 12,0 = ^ 22,0 = 0. Present nonequilibrium solution. ; 

differential Reynolds stress equation, . 

strain parameter is assumed to vary slowly, and an explicit expression can be found for the 
time evolution of the anisotropy of the Reynolds stress tensor in all planar homogeneous flow. 
The present nonequilibrium solution is extremely effective at capturing the initial behavior 
of the modeled Reynolds stress evolution, as well as the equilibrium states. In most cases, 
the present explicit nonequilibrium solution predicts stress anisotropies that are quite close 
to those given by the modeled differential Reynolds stress anisotropy evolution equation 
for all times; the small differences are attributed to the assumption of slow variation of 
the relative strain parameter used in obtaining the explicit expression of the time evolution 
of the modeled stress anisotropies. It has also been shown that the present nonequilibriuni 
solution is able to predict all the dynamic features of the Reynolds stress evolution, including 
the oscillatory nature of the stress anisotropy for elliptic flows. 
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APPENDIX A: REPRESENTATION OF b 2 - ±{b 2 }I 


Consider a symmetric, traceless tensor b for which the elements in any rectangular 
coordinate system are functions of the elements of two independent traceless tensors S* 
(symmetric) and W* (antisymmetric) in the same coordinate system, which is written as 

bij = k^Wu). 

The forms of these functional relationships also must be independent of the particular co- 
ordinate system in which they are expressed; that is, the relation between b, S*, and W* is 
isotropic. 8 

For two-dimensional mean flows, S* has one vanishing eigenvalue, and in the principal 
coordinate system of S*. the vort icily vector is aligned with the eigenvector of S* that 
corresponds to the vanishing eigenvalue. If the tensor b* is also assumed to have one 
eigenvector aligned with the eigenvector of S* that corresponds to t he vanishing eigenvalue, 
then the tensor b can be represented in terms of the tensors S* and W* and the scalar 
invariants {bS*}, {bW*S*}, and {bS* 2 }, as 5 

b= {bS*}S* + { b ^ -. S -l h s»W* - W*s*) + (i{bS* 2 }(S* 2 - fl). (Al) 

The quadratic term b 2 — ^ { b 2 } I can also be represented in terms of t he tensors S* and W~ 
and the scalar invariants {bS*}, {bW*S*}, and {bS* 2 }. 

If in expression (Al) the symmetric, traceless tensor b is replaced by b 2 — ^ {b 2 }I, then 
the following equation is obtained: 

I fb 2 W*S*f 

b 2 - -{b 2 }I = {b 2 S*}S* + 1 ~ ( S * W * — w * s *) (A *2) 

{W* z } 

+ f>({b 2 S- 2 )-i)b 2 ))(S- 2 -jI). 

Now, the scalar invariants {b 2 S*}, {b 2 W*S*}, and {b 2 S* 2 } in ( A2) must be expressed in 
terms of the scalar invariants {bS*}, {bW'S"}, and {bS* 2 }. 

For conciseness, relation (Al) can be rewritten as 

b = S«.'T il (AS) 


where the scalar coefficients a t are 

a i = {bS*}, «2 = {bW’S* } / { W* 2 }, « 3 = 6{bS* 2 }, 


and the tensors T, are given by 

T, = S\ T 2 = S*W* - W*S*, T : , = S’ 2 - -l. 


Therefore. 



{b 2 S*} = {b 2 T,}, {b 2 W*S*} = -|{b 2 T 2 }, {b 2 S* 2 } - ^{b 2 } = {b 2 T 3 }, 

and if (A3) is inserted into the above expressions, then 

{b 2 T i} = EE a i fl fc{ T i T * T *}’ (* = 1,2,3). (A4) 

j=l A;=l 


Finally, the 27 invariants {TjTyT;}, = 1,2,3) must be evaluated. As a result of 

symmetry properties ((j, k, i) = (i,j,k) = (fc, *,j)), only 11 invariants must be computed: 
= (1,1,1), (1,1,2), (1.1,3), (1,2,2), (1,2,3), (1,3,2), (1,3,3), (2,2,2), (2,2,3), 
(2,3,3), and (3,3,3). With the generalized version of the Cayley- Hamilton theorem, 9 the 
only resulting nonzero invariants are 


{T 2 T 3 ) 


1 

6 ’ 


{T 2 T 3 } = -^{W* 2 }. 


^ = -W 


together with the invariants that result from the cyclic permutations of the indices. There- 
fore, the relations 


i«,« 3 = 2{bS*}{bS* 2 }, 


{b 2 S*} 

{b 2 W*S*} = ^« 2 « 3 {W- 2 } = 2{bW*S~}{bS“ 2 }, 


3 


1 


.2 ^ ..2 r 1 X 7 * 2 1 ^ ,.2 ^ 


{b 2 S* 2 } - -{b 2 } = - -« 2 {W* 2 } - -a 2 = ^{bS*} i 


1 {bW*S*} 2 
3 {W* 2 } 


(A5) 

{bS* 2 } 2 , 


lead to the desired expression of the quadratic term in (18). 
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APPENDIX B: DERIVATION OF {bS*}, {bW*S*}, AND {bS* 2 } EQUATIONS 

Starting from the tensor evolution equation for the Reynolds stress anisotropy (8), a 
system of three scalar ordinary differential equations in the three scalar unknowns { bS“ } , 
{bW*S*}, and {bS* 2 } can be derived. 

By multiplying relation (8) by S*, taking the trace of the equation, and using the results 
of Appendix A to express {b 2 S*} in terms of {bS*} and {bS* 2 }, the following equation is 
obtained: 

A*} = {bS-}-2« 3 {bS* 2 } +2« 2 {bW*S*} (B1 ) 

of* gij 

+— {bS*}{bS* 2 } - a, - {L’S*}. 

V 

Similarly, multiplying equation (8) by either W*S* or S* 2 and taking the trace ol the equa- 
tion leads to the following equations, respectively: 

{^-W*S*} = - — {bw*s*} -rtXfbS*} + — {bW*S*}{bS* 2 } - {L*W’*S*}, (152) 

at* gij // 2 // 

and 

{^S 2 -] = - — jbS' 2 } - i«,{bS-} + ^l{bS-) 2 (IW) 

at* grj A (> 

+ -^{bW*S*) 2 - — {bS* 2 } 2 - {L*S 2 *}. 

:5// C v 

In obtaining these two equations, the following relations are used: 

{bS* :i } = ^{bS*}. 2{bW* 2 S*} + {bW*S*W*} = -t^{bS'}, 

which are consequences of the ( •ayley-Hamilton theorem. 9 Because the velocity gradients 
have been assumed independent of time, the following can easily be verified: 

An = -^-{bs*}. 

V//* ; dt* x ; 

f^WS-) = ^(bWS-J, (1*4 1 

{4^s 2 *} = 2-{bs- 2 ). 

x dt* 1 (It* 1 

Equations (Bl), (B2), and (B8) lead, therefore, to the desired system of scalar ordinary 
differential equations for the invariants {bS*}. {bW*S*}, and {bS* 2 }, 
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-^-{bS*} = - — {bS*} - 2a 3 {bS* 2 } + 2a 2 {bW*S*} 
df gi] 

—ax — {L*S*} 4- — {bS*}{bS‘ 2 }, 

V 

^-{bW*S*} = - — {bW*S*} - «X{bS*} - {L*W*S*} (B5) 

dt* gt] rj z 

+ — {bW*S‘}{bS* 2 }, 

V 

-^{bS- 2 | = --{bS' 2 ) - ia 3 {bS-} - {L"S* 2 } + ^{bS-) 2 
at grj 3 6 7] 

+ — ?J{bW*S*} 2 - — {bS* 2 } 2 . 

07/ C , 2 1 ] 
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APPENDIX C: SOLUTION OF ANISOTROPY EVOLUTION EQUATION 


The change of variables 


B x = 

r , «3 * o 2 'R 2 

B 2 = — C, 77 

Cl 2 H 

n - <■ 

transforms system (21) into the quadratic system of ordinary differential equations 

C = (■>«,/. - £) c, 

o = ^2o (/' — — ^ o T // (. ■. 


• / /c 

t/’ = 2 at 


(CM) 

( 02 ) 

(C3) 


y ■ 0 — c.i\ 


(04) 

(cr>) 

(CO) 


where 


H - -a 2 - 2 d\Tl 


2'T>2 


System ((.14) ((1(5) is subjected to the following initial conditions: 

</’o = V’(0) = By.o. 

0o — ^(0) = 2(i2 Him ~~ 2a^fi^ t o, 

Co = C(0) = (|fl2°3^2,0 — 2rt2^ 2 /^3.o) / H. 


(C7 


In this system, the evolution of the two variables (.:> and <6 is independent of the evolution 
of the variable C- Therefore, the quadratic system of ordinary differential equations 


d) — — —0 + II t + 'iftco. 

V 

: li . 2 

it = — ( ' + 0 + 2o y — a I , 

»/ 

can be solved, and the evolution of C is deduced by integrating (C l ). 


(08) 

(C9) 


C(H = Co exp 


f 


2 m/'(s) - 


a 


>i(- 


ds 


(CIO) 


By integrating (C8), the evolution of o can be given as a function of t/\ 

<p(D - Hcfo f‘ 0( r ) c -/o r + ^f* [2«y(.q-r/'/(dRy ((qi) 

Jo 
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Bv introducing the transformation 


4 in = - 


i 

2 a 




(C12) 


equation (Cll) can now be rewritten as 


* (n - + 


u)(0)(H + 2a(j) 0 ) + (3H J -d-^ds 


(C13) 


and equation (C9) can be written in terms of u> as the following integro-differential equation: 

£ = ()- + (H + 2a a, - - 3H f ^4rl.s - w (0 )(H + 2atf> 0 ), (C14) 

7 / Tj i JO 7}(SJ 

where u;(0) can take any nonzero value and u-’(O) = — u;(0)(2aV , o — ftfi) o) = — u>(0)(2cd3i,o — 
(j /r/o), with r/o = 7/(0) as the initial value of the relative strain parameter. Equation (C14) 
is integrated with the following transformation 

= j ^hh 

^2 = ^1 =w/l/, 

4h 3 = ^2 = do/ri + U?(l/r/), 

4^ = 4> 3 = Cjjt] 4- 2tj(l /rj) + uj{1/t]). 


The functions 4>, are the solution of 
’I'l = 4*2, 

^ = ^ f ... . (015) 

4> 3 = + \ H + 2 a«! - - 4> 2 + - (0 - 27/) 4*3 - ^V 2 (0 )(H + H"), 

v V vj v v 

where the initial conditions are given by 4 , i(0) = 0, 4*2(0) / 0, and 4* 3 (0) = — 4* 2 (0)(2af?i,o— 
ft / 7/ () - 7/0/7/0); and H° = 2ad>o = 4a(a 2 #2,o — a.3#3,o)- 

In the general case for which the relative strain parameter varies in time, system ( C 1 5 ) 
is a linear system of ordinary differential equations with variable coefficients. In the case of 
slow variations of 7 /, the approximations 

2-0, -~0 ( C 1 6 ) 

7 / 7 / 

are valid, and the solution of the system of ordinary differential equations (Cl 5) yields 

4b (n = A'E zvf f Ar< ‘ + HH + H°% 

r= 1 

4» 2 (C) = A'£> r e' vt \ (Cl 7) 

r = l 

* 3 (0 = KJ2drKt Xrt \ 

r~ 1 
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where K = V J , 2(0)[(A 2 — A1HA3 — Ai)(Aa — A 2 )] 1 and 


f'r = U 2 - 2aBifl\ r ~(H + H°) } (X p - A g ), 


A 3 — A 2 Ai - A3 A2 - Ai 

Ai A 2 A 3 


(C19) 


In (CIS), the indices p and q are chosen such that c rpq — —1. The A r are eigenvalues that 
are obtained as roots of the following third-order characteristic polynomial: 


A 3 - - A 2 - (// + 2aai)A + -H = 0. 


(C20) 


Finally, in terms of the original variables B, (I") the explicit solution is 


_L t _ 'M£l ] 

2a //(/*) ’ 


a >'R 2 

W) = 


_ (1 M!) _ ^2(0) 1 » 2 (0) g 

'//(/*) ^2 (t*) ® 2 (/*)J ^ 


ft jM£) 

'/(f*) 'hU*) 


jMO) 

* 2 (n 


In ((’21), the initial condition > ( 0 ) 0 is arbitrary, and its value can, therefore, be taken 

as vJ/ 2 (0) = 1. 
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